Modeling immune-parasite interactions

Model of Th1-Th2-P interactions

Let’s consider a model with a dynamic parasite population that both grows logistically and suffers mortality (or loss of biomass) due to immune response (but only one cell-type).

\[\begin{align} \frac{dT_i}{dt} &= b_i + c_i P + \frac{s_i T_i^2}{S_i^2 + T_i^2}\frac{I_{ij}}{I_{ij}+T_j} - m T_i, \\ \frac{dP}{dt} &= b_P P \left(1-\frac{P}{K_P}\right) - a T_2 \end{align}\]

## -------------------------------------------------------------------------------
## phaseR: Phase plane analysis of one- and two-dimensional autonomous ODE systems
## -------------------------------------------------------------------------------
## 
## v.2.1: For an overview of the package's functionality enter: ?phaseR
## 
## For news on the latest updates enter: news(package = "phaseR")

However, there is an issue with the above model. To see what it is, note the dynamics of the parasite again: \[\begin{equation} \frac{dP}{dt} = b_P P \left(1-\frac{P}{K_P}\right) - a T_2 \end{equation}\] I have used a formulation here where the mortality rate of the parasite is actually independent of the parasite density. The implication of this is that, if the parasite is driven extinct by the immune response, the growth rate is not zero, but rather is negative, leading the system to blow up in an unfortunate way. You can see that in a simulation where I make the immune induction by the parasite imbalanced, so that the parasite provokes a stronger Th2 immune response than Th1 (Fig. @ref(fig:imbalancedImmune)).

## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.56561e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.56561e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.56561e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.28127e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.28127e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.28127e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.04857e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.04857e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 1.04857e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 17.5505, R2 = 8.58131e-16
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 17.5505
##

## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps

## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go

If the parasite’s mortality rate is -a*T_2, and the parasite provokes a stronger Th2 immune response, it will be driven extinct and the parasite growth rate will continue to be negative, leading the system to take on negative abundances of parasites and immune cells.

We need to modify the parasite dynamics model to the following: \[\begin{equation} \frac{dP}{dt} = b_P P \left(1-\frac{P}{K_P}\right) - a T_2 P \end{equation}\]

To get some sense of this model’s potential dynamics, it’s actually quite useful to do some additional playing around with isocline analyses. However, in this case, I’m going to look at the isoclines for \(P\) and \(T_2\), treating \(T_1\) as a constant. This gives some sense of what’s possible, in terms of the outcomes of the relevant immunoparasite dynamics. One \(P\) isocline, given by \(dP/dt=0\), is the line \[\begin{equation} T_2 = \frac{b_P (K_P-P)}{a K_P}, \end{equation}\] with intercepts at \((T_2=b_P/a, P=0)\) and \((T_2=0, P=K_P)\). The other is the line \(P=0\). The \(T_2\) isocline, given by \(dT_2/dt = 0\), is a cubic function. You can see a sample configuration for a balanced immune response in Fig. @ref(fig:T2Pisoclines). Here there are four intersections between the three nullclines, corresponding to two stable equilibria (a high parasite, low Th2 equilibrium and a low parasite, high Th2 equilibrium), and two unstable equilibria (the equilibrium corresponding to zero parasites and “naive” immune response and an intermediate equilibrium).

Isoclines for T2 and P, assuming T1=200.

From this, however, it’s obvious that there are other possible ways that the isoclines can intersect. In particular, we can consider two ways of modifying the \(P\) nullcline by changing the parameters of the parasite growth equation that change the intersections. If I increase the carrying capacity for the parasite, then there are two equilibria, only one of which is stable, and represents the coexistence of both Th2 and P at a relatively high point (Fig. @ref(fig:T2Pisocline2).

Nullcline configuration and stability if the parasite carrying capacity is increased.

On the other hand, if increase the immune system’s killing ability by increasing \(a\), then there are two equilibria, one of which is stable, and represents the coexistence of both Th2 and P at a low immune biomass (Fig. @ref(fig:T2Pisocline3).

Nullcline configuration and stability if the immune attack rate is increased.

The other way that the isocline configuration can be modified is by changing the shape of the Th2 isocline. In particular, if I lower the “trough” of the function so that it crosses the Th2 axis, there are other dynamical possibilities that appear. For example, if I reduce the immune loss rate \(m\) from 1 to 0.95, I can get this configuration, with 4 equilibria, one of which is stable and represents the stable coexistence of the immune system and parasite with a high Th2 biomass and a low parasite biomass (Fig. @ref(fig:T2Pisocline4))

Nullcline configuration and stability when immune cell loss rate is lower.

If I take that configuration but reduce the carrying capacity of the parasite, then I get the configuration in Fig. @ref(fig:T2Pisocline5), with six equilibria, two of which are stable and represent the coexistence of the immune system and parasite at either an immune-dominated equilibrium, or a parasite-dominated one.

Nullcline configuration and stability when immune cell loss rate is lower and the parasite carrying capacity is decreased.

If I take that configuration, but increase the immune system’s killing ability, then I get the configuration in Fig. @ref(fig:T2Pisocline6), with five equilibria, two of which are stable and represent the coexistence of the immune system and parasite at an immune-dominated equilibrium or a parasite extinction equilibrium. This is a particularly interesting case, becuase it represents the two most biologically plausible outcomes of infection.

Nullcline configuration and stability when immune cell loss rate is lower and the parasite carrying capacity is decreased, and the immune killing rate is increased.

Finally, if I increase the parasite’s carrying capacity, I get the configuration in Fig. @ref(fig:T2Pisocline7), with three equilibria, only one of which is plausible and represents parasite extinction. Contrasting this with the results in Fig. @ref(fig:T2Pisocline6) is quite interesting, because it indicates that increasing the parasite’s carrying capacity can actually put the system into a state where the only possible outcome is an acute infection. If we think about the model as a model of parasite biomass accumulation, then the total biomass that is possible should actually depend on dose in an important way. In other way, the maximum parasite biomass should increase with dose as it increases the number of worms that could potentially reach adulthood. Thus, increasing dose (and increasing parasite biomass carrying capacity) can shift the configuration of the system from one where the parasite can persist to one where it cannot.

Nullcline configuration and stability when immune cell loss rate is lower and the parasite carrying capacity is decreased, and the immune killing rate is increased.

One way to build an extension to the model would be to more carefully consider the two aspects of parasite growth within the host, which is in terms of the numbers of parasites and the biomass of those parasites. The immune system can act to reduce the number of parasites, and it can act to decrease the biomass growth rate of the parasites. This may be an important extension to consider that more carefully deals with dose in a mechanistic way.

However, it is important to note that in all of the above, I was only considering a single Th1 abundance (\(T_1=200\)). What happens if the parasite is also secreting Th1 cytokines? How does the Th2 and P nullcline configuration change as the T1 cell population grows?

I will take the configuration shown in Fig. @ref(fig:T2Pisocline6) and increase the Th1 abundance. You can see that, as the Th1 cell population grows, eventually the acute infection equilibrium vanishes. This suggests that there may be two ways for the parasite to succeed. One is by inducing a strong enough Th1 response to eliminate the possibility of an acute infection. The other is to avoid inducing a strong immune response at all.

The challenge here is figuring out what to plot. What things are likely to vary among host or parasite genotypes? One thing that’s fairly easy to do is to look at the parasite isocline to clarify how changing those parameters will affect things. Recall that the intercepts of this linear equation are \((T_2=b_P/a, P=0)\) and \((T_2=0, P=K_P)\). Thus, changing the carrying capacity for the parasite slides the isocline up the y-axis, keeping the x-intercept fixed. The effects of this change will depend, of course, on the relative position of the \(T_2\) isocline. If the \(T_2\) isocline intercepts the \(T_2\) axis only once (near 0), then increasing the carrying capacity will tend to stabilize high immune, high parasite equilibrium; reducing the carrying capacity will tend to stabilize the low immune, low parasite equilibrium; intermediate carrying capacities can be bistable between these two equilibria. If the \(T_2\) isocline intercepts the \(T_2\) axis three times, the the effect of changing the carrying capacity depends quite a lot on where the \(P\) isocline intercepts the \(T_2\) axis. If the \(P\) isocline intercepts the \(T_2\) axis beyond the third intercept of the \(T_2\) isocline, then increasing the carrying capacity will tend to cause the system to stabilize at a high \(T_2\), low parasite biomass equilibrium. Lowering the carrying capacity will cause the system to become bistable between a high \(T_2\) and a low \(T_2\) equilibrium. If the \(P\) isocline intercepts the \(T_2\) axis between the second and third intercepts of the \(T_2\) isocline, then the acute infection equilibrium will always be stable, and increasing the carrying capacity will cause that to be the only stable equilibrium, whereas lowering the carrying capacity will cause bistability. If the \(P\) isocline intercepts the \(T_2\) axis between the first and second intercepts of the \(T_2\) isocline, then there will always be bistability between an acute and chronic infection and increasing the carrying capacity will just increase the \(P\) and \(T_2\) abundances at the chronic equilibrium.

Of course, it is reasonable to argue that carrying capacity might be fairly irrelevant to most infections, which will never get anywhere near their carrying capacities any way. In such a case it is more reasonable to ask what might happen as you change either the parasite’s growth rate \(b_P\) or the immune response’s attack rate \(a\), which together determine the \(T_2\) intercept (\(b_P/a\)). Again, this will depend on the number of intercepts of the \(T_2\) isocline with the \(T_2\) axis. If there is a single intercept, then increasing \(b_P/a\) will tend to cause the system to become stable at a high \(P\), high \(T_2\) equilibrium, whereas decreasing it will cause the system to become stable at a low \(P\), low \(T_2\) equilibrium. If there are three intercepts, then increasing \(b_P/a\) will have similar effects. However, if \(b_P/a\) is small enough (lower than the third intercept of the \(T_2\) isocline), then the acute infection equilibrium will always be stable, a key point.

Thus, understanding the qualitative dynamical possibilities we are interested in (acute versus chronic) requires answering two questions: (1) how many times does the \(T_2\) isocline intercept the \(T_2\) axis; (2) if that answer is 3, is \(b_P/a\) smaller than the third intercept? Unfortunately, there is essentially no analytical tractability to that question, as we are dealing with a cubic function. The zeros of a cubic function are complicated and, with a model of this complexity, there are no nice simplifications that make it possible to figure out how many there will be, unfortunately.

The best thing to do in this scenario is just play out how changing the various parameters changes the shape of the isocline.

Increasing the background immune loss rate increases the height of the peak and the trough of the cubic function. This decreases the chances of a stable acute infection equilibrium and makes it more likely that the only stable equilibrium will be one where the parasite persists at a high abundance, with a low abundance of Th2 cells.

Increasing the parasite’s stimulation of \(T_2\) decreases the heights of the peak and trough, but without affecting the intercepts of the isocline with the \(T_2\) axis. Thus, if the acute infection equilibrium is stable (because of the position of the \(P\) isocline), increasing the stimulation of the \(T_2\) response will tend to destabilize the chronic equilibrium, as you would expect. However, if there is bistability between two chronic equilibria (a high parasite and low parasite equilibrium), then you will end up in a situation where only the high \(T_2\), low \(P\) equilibrium is stable. Again, this is exactly what you would expect.

Increasing the autostimulation of the \(T_2\) response by increasing \(s_2\) drastically lowers the peak and trough of the \(T_2\) isocline. When \(s_2\) is small, the only stable equilibrium is one where \(P\) is relatively high and \(T_2\) is low. When \(s_2\) is large, the only stable equilibrium is the acute infection equilibrium. Again, this is exactly what you would expect.

Increasing the cross-inhibition of \(T_2\) by \(T_1\) has a relatively small effect (potentially because of the choice of \(T_1\)), and this effect is primarily on lowering the trough of the \(T_2\) isocline. Thus increasing the cross-inhibition could actually make an acute infection more likely (e.g., moving from bistability between two chronic equilibria bistability between an acute an a chronic equilibrium). That is the only response that’s a bit unexpected, although it does explain why bistability requires a large value of the cross-inhibition parameter.

I have struggled to find situations where increasing dose can move an infection from acute to chronic. The opposite is fairly easy to find. I think this is because the model I was using was not incorporating the lessons I learned from the analysis of simpler models. In particular, I have induction of immunity as a linear function of parasite biomass, whereas the previous simple model showed that bistability required a nonlinear (saturating) function. Here is the new model: \[\begin{align} \frac{dT_i}{dt} &= b_i + c_i \frac{P}{C_i + P} + \frac{s_i T_i^2}{S_i^2 + T_i^2}\frac{I_{ij}}{I_{ij}+T_j} - m T_i, \\ \frac{dP}{dt} &= b_P P \left(1-\frac{P}{K_P}\right) - a T_2 P \end{align}\]

You can see that bistability in the Th1-Th2 plane is much more impervious to changes in parasite biomass in this model.

bistability in the Th1-Th2 plane is much more impervious to changes in parasite biomass in this model.

Here you can see how changing dose alters the configuration of the isoclines, showing the Th1 and Th2 isoclines for a low dose (black and red) and a high dose (blue and green) and initial Th1-Th2 densities where a low dose leads to a chronic infection but a high dose leads to clearance. The flows show that, from an initial low dose, the movement in the phase plane is towards increasing Th1-ness and decreasing Th2-ness (grey) whereas an initial high dose flows towards increasing Th1-ness and increasing Th2-ness. The full picture is shown by pairing the Th1-Th2 isoclines with the Th2-P isoclines. You can see that, from a low dose, the parasite biomass is initially increasing, whereas from a high dose, the parasite biomass is initially decreasing. Putting the pictures together: from a low dose, the parasite biomass is increasing and the Th1 response is increasing while the Th2 response is decreasing (due to the fact that the system starts with a Th1 bias that allows the Th1 response to suppress the Th2); at a high dose, the parasite bimoass is decreasing and the Th1 and Th2 response are both increasing (due to the fact that the parasite population is larger, which stimulates the Th2 response).

Here are the isoclines for a low dose (black and red) and a high dose (blue and green) and initial Th1-Th2 densities where a low dose leads to a chronic infection but a high dose leads to clearance. The flows show that, from an initial low dose, the movement in the phase plane is towards increasing Th1-ness and decreasing Th2-ness (grey) whereas an initial high dose flows towards increasing Th1-ness and increasing Th2-ness.

Figure for grant proposal.

Figure showing how initial Th2-ness (x-axis) influences when infections are chronic (left of the solid line) versus chronic (right of the solid line), for different initial parasite doses (y-axis) and different initial numbers of T cells (panels). In this figure, the rate of Th1 induction by the parasite is 200 and the rate of Th2 induction is 0.